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Abstract — This paper presents control and coordination al- 
gorithms for groups of vehicles. The focus is on autonomous 
vehicle networks performing distributed sensing tasks where 
each vehicle plays the role of a mobile tunable sensor. The 
paper proposes gradient descent algorithms for a class of 
utility functions which encode optimal coverage and sens- 
ing policies. The resulting closed-loop behavior is adaptive, 
distributed, asynchronous, and veriflably correct. 

Keywords — Coverage control, distributed and asyn- 
chronous algorithms, centroidal Voronoi partitions 



I. Introduction 
Mobile sensing networks 

The deployment of large groups of autonomous vehi- 
cles is rapidly becoming possible because of technological 
advances in networking and in miniaturization of electro- 
mechanical systems. In the near future large numbers of 
robots will coordinate their actions through ad-hoc com- 
munication networks and will perform challenging tasks 
including search and recovery operations, manipulation in 
hazardous environments, exploration, surveillance, and en- 
vironmental monitoring for pollution detection and esti- 
mation. The potential advantages of employing teams of 
agents are numerous. For instance, certain tasks are diffi- 
cult, if not impossible, when performed by a single vehicle 
agent. Further, a group of vehicles inherently provides ro- 
bustness to failures of single agents or communication links. 

Working prototypes of active sensing networks have al- 
ready been developed; sec [Q, ||. In ||, launchable 
miniature mobile robots communicate through a wireless 
network. The vehicles are equipped with sensors for vibra- 
tions, acoustic, magnetic, and IR signals as well as an active 
video module (i.e., the camera or micro-radar is controlled 
via a pan-tilt unit). A second system is suggested in j| 
under the name of Autonomous Oceanographic Sampling 
Network; see also §, @. In this case, underwater vehi- 
cles are envisioned measuring temperature, currents, and 
other distributed oceanographic signals. The vehicles com- 
municate via an acoustic local area network and coordinate 
their motion in response to local sensing information and to 
evolving global data. This mobile sensing network is meant 
to provide the ability to sample the environment adaptively 
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in space and time. By identifying evolving temperature 
and current gradients with higher accuracy and resolution 
than current static sensors, this technology could lead to 
the development and validation of improved oceanographic 
models. 

Optimal sensor allocation and coverage problems 

A fundamental prototype problem in this paper is that of 
characterizing and optimizing notions of quality-of-service 
provided by an adaptive sensor network in a dynamic en- 
vironment. To this goal, we introduce a notion of sensor 
coverage that formalizes an optimal sensor placement prob- 
lem. This spatial resource allocation problem is the sub- 
ject of a discipline called locational optimization 0, ||, 

Locational optimization problems pervade a broad spec- 
trum of scientific disciplines. Biologists rely on locational 
optimization tools to study how animals share territory and 
to characterize the behavior of animal groups obeying the 
following interaction rule: each animal establishes a region 
of dominance and moves toward its center. Locational opti- 
mization problems arc spatial resource allocation problems 
(where to place mailboxes in a city or cache servers on the 
internet) and play a central role in quantization and infor- 
mation theory (the design of a minimum-distortion fixed- 
rate vector quantizer is a locational problem) . Other tech- 
nologies affected by locational optimization include mesh 
and grid optimization methods, clustering analysis, data 
compression, and statistical pattern recognition. 

Because locational optimization problems are so widely 
studied, it is not surprising that methods are indeed avail- 
able to tackle coverage problems; see 0, |Q, [Q, [O. 
However, most currently-available algorithms are not ap- 
plicable to mobile sensing networks because they inherently 
assume a centralized computation for a limited size prob- 
lem in a known static environment. This is not the case in 
multi-vehicle networks which, instead, rely on a distributed 
communication and computation architecture. Although 
an ad-hoc wireless network provides the ability to share 
some information, no global omniscient leader might be 
present to coordinate the group. The inherent spatially- 
distributed nature and limited communication capabilities 
of a mobile network invalidate classic approaches to algo- 
rithm design. 

Distributed asynchronous algorithms for coverage control 

In this paper we design coordination algorithms imple- 
mentable by a multi-vehicle network with limited sensing 
and communication capabilities. Our approach is related 
to the classic Lloyd algorithm from quantization theory; 
see [^3| for a reprint of the original report and juj for a 
historical overview. We present Lloyd descent algorithms 
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that take into careful consideration all constraints on the 
mobile sensing network. In particular, we design coverage 
algorithms that are adaptive, distributed, asynchronous, 
and verifiably asymptotically correct: 

Adaptive: Our coverage algorithms provide the network 
with the ability to address changing environments, sens- 
ing task, and network topology (due to agents departures, 
arrivals, or failures). 

Distributed: Our coverage algorithms arc distributed in the 
sense that the behavior of each vehicle depends only on 
the location of its neighbors. Also, our algorithms do 
not required a fixed-topology communication graph, i.e., 
the neighborhood relationships do change as the network 
evolves. The advantages of distributed algorithms are scal- 
ability and robustness. 

Asynchronous: Our coverage algorithms are amenable to 
asynchronous implementation. This means that the al- 
gorithms can be implemented in a network composed of 
agents evolving at different speeds, with different compu- 
tation and communication capabilities. Furthermore, our 
algorithms do not require a global synchronization and con- 
vergence properties are preserved even if information about 
neighboring vehicles propagates with some delay. An ad- 
vantage of asynchronism is a minimized communication 
overhead. 

Verifiable Asymptotically Correct: Our algorithms guaran- 
tees monotonic descent of the cost function encoding the 
sensing task. Asymptotically the evolution of the mobile 
sensing network is guaranteed to converge to so-called cen- 
troidal Voronoi configurations that are critical points of the 
optimal sensor coverage problem. 

Let us describe in some detail what are the contribu- 
tions of this paper. Section || reviews certain locational 
optimization problems and their solutions as centroidal 
Voronoi partitions. Section provides a continuous-time 
version of the classic Lloyd algorithm from vector quantiza- 
tion and applies it to the setting of multi- vehicle networks. 
In discrete-time, we propose a family of Lloyd algorithms. 
We carefully characterize convergence properties for both 
continuous and discrete-time versions (Appendix VII col- 
lects some relevant facts on descent flows). We discuss a 
worst-case optimization problem, we investigate a simple 
uniform planar setting, and we present numerical results. 



Section IV presents two asynchronous distributed imple- 



mentations of Lloyd algorithm for ad-hoc networks with 
communication and sensing capabilities. Our treatment 
carefully accounts for the constraints imposed by the dis- 
tributed nature of the vehicle network. We present two 
asynchronous implementations, one based on classic results 
on distributed gradient flows, the other based on the struc- 
ture of the coverage problem. 



Section V-A considers vehicle models with more realistic 
dynamics. We present two formal results on passive ve- 
hicle dynamics and on vehicles equipped with individual 
local controllers. We present numerical simulations of pas- 
sive vehicle models and of unicycle mobile vehicles. Next, 



Section V-B describes density functions that lead the multi- 



vehicle network to predetermined geometric patterns. 



Review of distributed algorithms for cooperative control 

Recent years have witnessed a large research effort fo- 
cused on motion planning and coordination problems for 
multi- vehicle systems. Issues include geometric patterns 
U, 01, formation control @, |§, |§, |o), §l|, and 
conflict avoidance p2| , |2j| . Algorithms for robotic sens- 
ing tasks are presented for example in |24|] , p5fl . It is only 
recently, however, that truly distributed coordination laws 
for dynamic networks are being proposed; ej?-, see pq ], p7[ 
and the conference versions of this work P8[ , |29f |. 

Heuristic approaches to the design of interaction rules 
and emerging behaviors have been throughly investigated 
within the literature on behavior-based robotics; see j3(J, 
H, 0, |1, @ Jp, @. An example of coverage 
control is discussed in [|37|. Along this line of research, al- 
gorithms have been designed for sophisticated cooperative 
tasks. However, no formal results are currently available 
on how to design reactive control laws, ensure their cor- 
rectness, and guarantee their optimality with respect to an 
aggregate objective. 

The study of distributed algorithms is concerned with 
providing mathematical models, devising precise specifica- 
tions for their behavior, and formally proving their cor- 
rectness and complexity. Via an automata-theoretic ap- 
proach, the references [Q, [[59) treat distributed consensus, 
resource allocation, communication, and data consistency 
problems. From a numerical optimization viewpoint, the 
works in |4(J , (4^| , |^] discuss distributed asynchronous al- 
gorithms as networking algorithms, rate and flow control, 
and gradient descent flows. Typically, both these sets of 
references consider networks with fixed topology, and do 
not address algorithms over ad-hoc dynamically changing 
networks. Another common assumption is that, any time 
an agent communicates its location, it broadcasts it to ev- 
ery other agent in the network. In our setting, this would 
require a non-distributed communication set-up. 

II. From location optimization to centroidal 
Voronoi partitions 

A. Locational optimization 

In this section we describe a collection of known facts 
about a meaningful optimization problem. References in- 
clude the theory and applications of centroidal Voronoi 
partitions, see [Q, and the discipline of facility location, 
sec [||. Along the paper, we interchangeably refer to the 
elements of the network as sensors, agents, vehicles, or 
robots. 

Let Q be a convex polytope in R N and let || ■ || denote the 
Euclidean distance function. We call a map <fi : Q — > R + 
a distribution density function if it represents a measure 
of information or probability that some event take place 
over Q. In equivalent words, we can consider Q to be the 
bounded support of the function 4>. Let P = (pi, ■ ■ ■ ,p n ) 
be the location of n sensors, each moving in the space Q. 
Because of noise and loss of resolution, the sensing per- 
formance at point q taken from ith sensor at the position 
Pi degrades with the distance \\q — pi\\ between q and pf, 
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we describe this degradation with a non-decreasing differ- 
cntiable function / : K + — ► R+. Accordingly, / (\\q — pi\\) 
provides a quantitative assessment of how poor the sensing 
performance is. 

Remark II. 1: As an example, consider n mobile robots 
equipped with microphones attempting to detect, identify, 
and localize a sound-source. How should we plan to robots ' 
motion in order to maximize the detection probability? As- 
suming the source emits a known signal, the optimal detec- 
tion algorithm is a matched filter (i.e., convolve the known 
waveform with the received signal and threshold). The 
source is detected depending on the signal-to-noisc-ratio, 
which is inversely proportional to the distance between the 
microphone and the source. Various electromagnetic and 
sound sensors have signal-to-noise ratios inversely propor- 
tional to distance. 

Within the context of this paper, a partition of Q is a 
collection of n polytopes W = {W\, . . . , W n } with disjoint 
interiors whose union is Q. We say that two partitions 
W and W are equal if Wi and W/ only differ by a set of 
0-measure zero, for all i G {1, . . . , n}. 

We consider the task of minimizing the locational opti- 
mization function 



and oo-norm over Q = R m , and Voronoi diagrams can be 
defined over Riemannian manifolds; see ]4^| . Some useful 
facts about the Euclidean setting are the following: if Q is 
a convex polytope in a TV-dimensional Euclidean space, the 
boundary of each Vj is the union of (N — l)-dimensional 
convex polytopes. 

In what follows, we shall write 



H V (P) = H(P,V(P)). 



Note that 



H V (P) = [ min f(\\q- Pi \\)dcj>(q), 
Jn ie{l,...,n} 



(2) 



E 



min 

n} 



that is, the locational optimization function can be inter- 
preted as an expected value composed with a min opera- 
tion. This is the usual way in which the problem is pre- 
sented in the facility location and operations research lit- 
erature §!,{§. Remarkably, one can show that 



(1) 



dH v 

dpi 



(P) 



&H 

dpi 



(P,V(P)) 



d_ 

dp, 



||) #(<?), 



(3) 



and deduce some smoothness properties of Tiy. Since 
the Voronoi partition V depends at least continuously on 
P = (pi, . . . ,p n ), the function Hy is at least continuously 
diffcrcntiablc. 

C. Centroidal Voronoi partitions 

Let us recall some basic quantities associated to a region 
V C M. N and a mass density function p. The (generalized) 
mass, ccntroid (or center of mass), and polar moment of 
inertia are defined as 

AI V = / p(q)dq, C v = — - [ qp{q)dq, 
Jv M v Jv 



where we assume that the zth sensor is responsible for mea- 
surements over its "dominance region" Wi. Note that the 
function Ti. is to be minimized with respect to both (1) 
the sensors location P, and (2) the assignment of the dom- 
inance regions W. This problem is referred to as a fa- 
cility location problem and in particular as a continuous 
p-median problem in Q . 

Remark II. 2: Note that if we interchange the positions of 
any two agents, along with their associated regions of dom- 
inance, the value of the locational optimization function Tt 
is not affected. To eliminate this discrete redundancy, one 
could take the discrete group of permutations E n with the 
natural action on Q n , and consider Q n /T, n as the configu- 
ration space for the position P of the n vehicles. 

B. Voronoi partitions 

One can easily see that, at fixed sensors location, the 
optimal partition of Q is the Voronoi partition V(P) = 
{Vi, . . . , V n ] generated by the points (pi, . . . ,p n ): 

Vi = {qeQ \ \\q- Pi \\ < \\q-pj\\, Vj^i}. 

We refer to for a comprehensive treatment on Voronoi 
diagrams, and briefly present some relevant concepts. The 
set of regions {Vi, . . . ,V n } is called the Voronoi diagram 
for the generators {pi, ■ ■ ■ ,p n }- When the two Voronoi 
regions Vi and Vj are adjacent, pi is called a (Voronoi) 
neighbor of pj (and vice- versa) . The set of indexes of the 
Voronoi neighbors of pi is denoted by N(i). Clearly, j e 
Af(i) if and only if i G Af(j). We also define the (i,j)-face 

as Aij = Vi n Vj. Voronoi diagrams can be defined with that is, we assume — Pi\\) = \\q — Pi\\ 2 - Applying the 
respect to various distance functions, e.g., the 1-, 2-, s-, parallel axis theorem leads to simplifications for both the 



J v,p = / \\q-p\\ p{q)dq. 

Jv 

Additionally, by the parallel axis theorem, one can write, 

Jv, P = J v> c v +M v \\p-C v \\ 2 (4) 

where Jv,c v € ^+ i s defined as the polar moment of inertia 
of the region V about its centroid. 

Let us consider again the locational optimization prob- 
lem (Q), and suppose now we are strictly interested in the 
setting 



n 

H(P,W)=y;/ \\q-Pi\\ 2 dct>(q), 

„-_i JWi 



(5) 
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function 7iy and its partial derivative: 



dUy , 

dpi 



i=l i=l 

-(P) = 2M Vi (pi-C Vi ). 



It is convenient to define H\> 1 = X)"=i *^Vj c v an d Tiy 2 = 

E^iMvJ|ft-c Vi || a . 

Therefore, the (not necessarily unique) local minimum 
points for the location optimization function 7iy arc cen- 
troids of their Voronoi cells, i.e., each location pi satisfies 
two properties simultaneously: it is the generator for the 
Voronoi cell V% and it is its centroid 

Cvt = argmin p . Hy(P). 

Accordingly, the critical partitions and points for H are 
called centroidal Voronoi partitions. We will refer to a sen- 
sors configuration as a centroidal Voronoi configuration if 
it gives rise to a centroidal Voronoi partition. This discus- 
sion provides a proof alternative to the one given in Jljj for 
the necessity of centroidal Voronoi partitions as solutions 
to the continuous p-mcdian location problem. 

iii. continuous and dlscrete-tlme lloyd 
Descent for Coverage Control 

In this section, we describe algorithms to compute the 
location of sensors that minimize the co st H, both in con- 
tinuous and in discrete-time. In Section [II- A, we propose 



a continuous-time version of the classic Lloyd algorithm. 
Here, both the positions and partitions evolve in continu- 
ous time, whereas Lloyd algorithm for vector quantization 
is designed in discrete time. In Section III-B, we develop 



a family of variations of Lloyd algorithm in discrete time. 
In both setting, we prove that the proposed algorithms are 
gradient descent flows. 

A. A continuous-time Lloyd algorithm 

Assume the sensors location obeys a first order dynami- 
cal behavior described by 



Pi 



Consider Tiy a cost function to be minimized and impose 
that the location pi follows a gradient descent. In equiv- 
alent control theoretical terms, consider TL\> a Lyapunov 
function and stabilize the multi-vehicle system to one of 
its local minima via dissipative control. Formally, we set 



>(Pi ~ CvJ, 



(6) 



Q. Assuming this set is finite, the sensors location con- 
verges to a centroidal Voronoi configuration. 
Proof: Under the control law (||), we have 

l — l 

n 

= -2k pmp M Vi \\Pi - C Vi || 2 = -2k pIop n v , 2 (P(t)). 

i=l 

By LaSalle's principle, the sensors location converges to the 
largest invariant set contained in 7iy 2 (0); which is precisely 
the set of centroidal Voronoi configurations. Since this set 
is clearly invariant for (^), we get the stated result. If 
7iy2(0) consists of a finite collection of points, then P(t) 



converges to one of them, see Corollary VII. 2 

Remark III. 2: If Hy 1 (0) is finite, and P{t) -> C, then 
a sufficient condition that guarantees exponential conver- 
gence is that the Hessian of H\> be positive definite at C. 
This property is known to be an open problem, see 
Note that this gradient descent is not guaranteed to find 
the global minimum. For example, in the vector quantiza- 
tion and signal processing literature |l4j ], it is known that 
for bimodal distribution density functions, the solution to 
the gradient flow reaches local minima where the number 
of generators allocated to the two region of maxima are not 
optimally partitioned. 

B. A family of discrete-time Lloyd algorithms 

Let us consider the following class of variations of Lloyd 
algorithm. Let T be a continuous mapping T : Q n — > Q n 
verifying the following two properties: 

(a) for alii e {1, ... ,n}, \\Ti{P)-C Vi{P) \\ < \\ Pi -C Vi(P) l 
where 7j denotes the ith component of T, 

(b) if P is not centroidal, then there exists a j such that 
\\T 3 (P) - C V](P} \\ < \\p 3 -C V]{P) \\. 

Property (a) guarantees that, if moving, the agents of the 
network do not increase their distance to its correspond- 
ing centroid. Property (b) ensures that at least one robot 
moves at each iteration and strictly approaches the centroid 
of its Voronoi region. Because of this property, the fixed 
points of T are the set of centroidal Voronoi configurations. 

Proposition III. 3 (Discrete-time Lloyd descent) Let Pq 
G Q n denote the initial sensors location. Then, the 
sequence {r™ l (Po)}m>i converges to the set of cen- 
troidal Voronoi configurations. If this set if finite, then 
{T m (P )} m >i converges to a centroidal Voronoi configura- 
tion. 

Proof: Consider H\> : Q n — ► as an objective 
function for the algorithm T. Note that 

H(P,V(P))<H(P,W), (7) 
Moreover, the parallel 



with strict inequality if V(P) 
axis theorem guarantees 



where k is a positive gain, and where we assume that the 
partition V(P) — {Vi, . . . , V n } is continuously updated. 

Proposition III.l (Continuous-time Lloyd descent) For the 
closed- loop system induced by equation (||) , the sensors lo- 
cation converges asymptotically to the set of critical points as long as — C\Yi II < \\Pi — C-Wi 
of 7iy, i.e., the set of centroidal Voronoi configurations on with strict inequality if for any i, \\p[ 



H{P',W) < H{P,W), 



(8) 



for all i G {1, . . . , n}, 

-cwj < \\Pi-CwA\- 
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In particular, Tt(Cyv, W) < H(P, W), with strict inequality 
if P ^ Cw, where Cyy denotes the set of centroids of the 
partition W. 
Now, we have 

H V (T(P)) = H(T(P),V(T(P))) < H(T(P), V(P)) , 

because of M). In addition, because of property (a) of T, 
inequality (g) yields 

H(T(P), V(P)) < H(P, V(P)) = «v(P) , 

and the inequality is strict if P is not centroidal by property 
(b) of T. Hence, TLy is a descent function for the algorithm 
T. The result now follows fro m the global convergence 



Theorem VIL3 and Proposition |VII.4| . ■ 
Remark III. J.: Lloyd algorithm in quantization the- 
ory |l4[ is usually presented as follows: given the loca- 
tion of n agents, pi, . . . ,p n , (i) construct the Voronoi par- 
tition corresponding to P = (pi, . . . ,p n ); (ii) compute the 
mass centroids of the Voronoi regions found in step (i). Set 
the new location of the agents to these centroids; and re- 
turn to step (i) . Lloyd algorithm can also be seen as a fixed 
point iteration. Consider the mappings LLi : Q n —> Q for 
ie{l,...,n} 



LLi(pi,. ..,p„) 



'Vi(P) 



(j)(q)dq 



q<p(q)dq . 



Vi{P) 



Let LL : Q n -> Q n be defined by LL = (LL U . . . , LL n ). 
Clearly, LL is continuous (indeed, C 1 ), and corresponds to 
Lloyd algorithm. Now, \\LL t (P) - C Vi || = < \\p { - Cy. ||, 
for all i G {1, ...,n}. Moreover, if P is not centroidal, 
then the inequality is strict for all pi ^ Cy ■ Therefore, LL 
verifies properties (a) and (b). 

C. Generalized settings, worst-case design, and the p- 
center problem 

Different sensor performance functions / in equation (^) 
correspond to different optimization problems. Provided 
one uses the Euclidean distance in the definition of TLw, 
the standard Voronoi partition computed with respect to 
the Euclidean metric remains the optimal partition. For 
arbitrary /, it is not possible anymore to decompose Tiy 
into the sum of terms similar to 7iy,i and 7iy.2- Neverthe- 
less, it is still possible to implement the gradient flow via 
the expression for the partial derivative (y). 

Proposition III. 5: Assume the sensors location obeys a 
first order dynamical behavior, pi = tt.j. Then, for the 
closed- loop system induced by the gradient law (||), Ui = 
—dHv/dpi, the sensors location P = (pi, . . . ,p n ) converges 
asymptotically to the set of critical points of 7iy . Assuming 
this set is finite, the sensors location converges to a critical 
point. 

More generally, various distance notions can be used to 
define locational optimization functions. Different perfor- 
mance function gives rise to corresponding notions of "cen- 
ter of a region" (any notion of geometric center, mean, 
or average is an interesting candidate). These can then 



be adopted in designing coverage algorithms. We refer 
to |44| l for a discussion on Voronoi partitions based on non- 
Euclidean distance functions and to , jl^] for a discussion 
on the corresponding locational optimization problems. 

Next, let us discuss an interesting variation of the origi- 
nal problem. In Q , Q , minimizing the expected minimum 
distance function 7Yy in equation (^) is referred to as the 
continuous p-median problem. It is instructive to consider 
the worst-case minimum distance function, corresponding 
to the scenario where no information is available on the 
distribution density function. In other words, the network 
seeks to minimize the largest possible distance from any 
point in Q to any of the sensor locations, i.e., to minimize 
the function 



max mm 



II? -Pill 



max 

i€{l,...,n} 



max o • 



Pi II 



This optimization is referred to as the p- center problem 
in H , . One can design a strategy for the p-center prob- 
lem analog to the Lloyd algorithm for the p-median prob- 
lem: each vehicle moves, in continuous or discrete-time, 
toward the center of the minimum-radius sphere enclosing 
the polytope. To the best of our knowledge, no conver- 
gence proof is available in the literature for this algorithm; 
e.g., see p5[ . We refer to for a convergence analysis of 
the continuous and discrete time algorithms. 

In what follows, we shall restrict our attention to the 
p-median problem and to centroidal Voronoi partitions. 



D. Computations over polygons with uniform density 

In this section, we investigate closed-form expression for 
the control laws introduced above. Assume the Voronoi 
region Vi is a convex polygon (i.e., a polytope in R 2 ) with 
Ni vcrtcxes labeled {(#0, Uo), ■■• , (xNi-x, UNi-i)} such as 
in Figure [|. It is convenient to define (a;^ , y^r. ) = (xo,yo). 
Furthermore, we assume that the density function is <p(q) = 
1. By evaluating the corresponding integrals, one can ob- 









(x2,y-i) r 




• / 






(x ,y ) = 








Fig 


1 



Notation conventions for a convex polygon. 
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tain the following closed-form expressions 

Ni-l 

2 



My i = O ^2 ( X kVk+l - Xk+lVk) 



k=0 



Ni-l 



Y~ ^2 ( Xk + x k+l)( x kVk+i - x k+ iy k ) (9) 



6M Vt 



k=0 
Ni-l 



Cv,, y = 53 (Vk + yk+i){x k yk+i 



Xk+iDk) 



k=0 



To present a simple formula for the polar moment of inertia, 
let x k = Xk-C Vi , x and y k = y k ~C Vl , y , for k € {0, . . . , JV»- 
1}. Then, the polar moment of inertia of a polygon about 
its centroid, Jvi,c becomes 



Ni-l 



J Vi,C Vi = Yl ( 



XkVk+i - Xk+iVk) 



k=0 

{x\ + XkXk+i + x 2 k+1 +yl + ykVk+i + J/fc+i) • 

The proof of these formulas is based on decomposing the 
polygon into the union of disjoint triangles. We refer to Q 
for analog expressions over R^. 

A second observation is that the Voronoi polygon's ver- 
texes can be expressed as a function of the neighboring 
vehicles. The vertexes of the ith Voronoi polygon which lie 
in the interior of Q are the circumcenters of the triangles 
formed by pi and any two neighbors adjacent to pi. The 
circumcenter of the triangle determined by Pi, pj, and pk 
is 



1 

4M 



\a k j\\ 2 (aji ■ a ik )pi + \\a lk \\ 2 (a kl ■ a^Pj 



+ llojill^aifc • a k j)pkj , (10) 



where M is the area of the triangle, and ai s = pi — p s - 

Equation (^|) for a polygon's centroid and equation ( [To| ) 
for the Voronoi cell's vertexes lead to a closed-form alge- 
braic expression for the control law in equation (||) as a 
function of the neighboring vehicles' location. 

E. Numerical simulations 

To illustrate the performance of the continuous-time 
Lloyd algorithm, we include some simulation results. The 
algorithm is implemented in Mathematica as a single cen- 
tralized program. For the R 2 setting, the code computes 
the bounded Voronoi diagram using the Mathematica pack- 
age ComputationalGeometry, and computes mass, cen- 
troid, and polar moment of inertia of polygons via the nu- 
merical integration routine NIntegrate. Careful attention 
was paid to numerical accuracy issues in the computation of 
the Voronoi diagram and in the integration. We illustrate 
the performance of the closed-loop system in Figure ||. 

IV. Asynchronous distributed implementations 

In this section we show how the Lloyd gradient algo- 
rithm can be implemented in an asynchronous distributed 



fashion. In Section [V-A we describe our model for a dis- 
tributed asynchronous network of robotic agents. Next, 
we provide two distributed algorithms for the local com- 
putation and maintenance of the Voronoi cells. Finally, in 
Section IV-C we propose two distributed asynchronous im- 



plementations of Lloyd algorithm: the first one is based on 
the gradient optimization algorithms as described in p0| 
and the second one relies on the special structure of the 
coverage problem. 

A. Modeling an asynchronous distributed network of mo- 
bile robotic agents 

We start by modeling a robotic agent that performs sens- 
ing, communication, computation, and control actions. We 
are interested in the behavior of the asynchronous net- 
work resulting from the interaction of finitely many robotic 
agents. A theoretical framework to formalize the following 
concepts is that developed in the theory of distributed al- 
gorithms; see H^]. 

Let us here introduce the notion of robotic agent with 
computation, communication, and control capabilities as 
the ith element of a network. The ith agent has a proces- 
sor with the ability of allocating continuous and discrete 
states and performing operations on them. Each vehicle 
has access to its unique identifier i. The ith agent occupies 
a location p, £ Q C and it is capable of moving in 
space, at any time t G M + for any period of time St E R+, 
according to a first order dynamics of the form: 



IMI<1, Vs€[t,t + dt}. 



The processor has access to the agent's location pi and 
determines the control pair (St,Ui). The processor of the 
ith agent has access to a local clock ti G M.+ U {0}, and 
a scheduling sequence, i.e., an increasing sequence of times 
{T hk G R+ U {0} k G N U {0}} such that T ifi = and 
U,min < T hk +i - Ti^ k < ^, max . The processor of the ith 
agent is capable of transmitting information to any other 
agent within a closed disk of radius Ri G M+ . We assume 
the communication radius Ri to be a quantity controllable 
by the ith processor and the corresponding communication 
bandwidth to be limited. 

We shall alternatively consider networks of robotic agents 
with computation, sensing, and control capabilities. In this 
case, the processor of the ith agent has the same compu- 
tation and control capabilities as before. Furthermore, we 
assume the processor can detect any other agent within a 
closed disk of radius Ri G R+. We assume the sensing 
radius Ri to be a quantity controllable by the processor. 

B. Voronoi cell computation and maintenance 

A key requirement of the Lloyd algorithms presented in 



Section III is that each agent must be able to compute 



its own Voronoi cell. To do so, each agent needs to know 
the relative location (distance and bearing) of each Voronoi 
neighbor. The ability of locating neighbors plays a central 
role in numerous algorithms for localization, media access, 
routing, and power control in ad-hoc wireless communica- 
tion networks; e.g., see Q, Q, Q, |3l) and references 
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Fig. 2 

Lloyd continuous-time algorithm for 32 agents on a convex polygonal environment, with Gaussian density 

<j> = exp(5.(— X 2 — y 2 )) CENTERED ABOUT THE GRAY POINT IN THE FIGURE. THE CONTROL GAIN IN (g) IS fc PRO p = 1 FOR ALL THE VEHICLES. 
THE LEFT (RESPECTIVELY, RIGHT) FIGURE ILLUSTRATES THE INITIAL (RESPECTIVELY, FINAL) LOCATIONS AND VORONOI PARTITION. THE 

CENTRAL FIGURE ILLUSTRATES THE GRADIENT DESCENT FLOW. 





Name: Adjust sensing radius algorithm 




Goal: distributed Voronoi cell 




Requires: sensor with radius Ri 


Local agent i performs: 


1 


initialize Ri, detect vehicles pj within radius Ri 


2 


update P l (ti), compute W (pi(ti) , Ri) 


3 


while Ri < 2max ge v^( Pi (t j ),fti) \\Pi(U) - q\\ do 


4 


Ri := 2max qeW ( Pi (t i ),R i ) \\Pi{U) - q\\ 


5 


detect vehicles Pj within radius Ri 


6 


update 


7 


compute W(pi(ti),Ri) 


8 


end while 


9 


set R4 := 2max,g ff ( ft ( t .)^j - q\\ 


10: set Vi := Wi(pi(ti),Ri) 



therein. Therefore, any motion control scheme might be 
able to obtain this information from the underlying com- 
munication layer. In what follows, we set out to provide 
a distributed asynchronous algorithm for the local compu- 
tation and maintenance of Voronoi cells. The algorithm is 
related to the synchronous scheme in |HJ and is based on 
basic properties of Voronoi diagrams. 

We present the algorithm for a robotic agent with sens- 
ing capabilities (as well as computation and control). The 
processor of the ith agent allocates the information it has 
on the position of the other agents in the state variable P\ 
The objective is to determine the smallest distance Ri for 
vehicle i which provides sufficient information to compute 
the Voronoi cell Vi. We start by noting that Vi is a subset 
of the convex set 



W(pi,Ri) = B( Pi ,Ri) n (n Mpi _ Pjll < Ri Sij), (li) 



where B(pi,Ri) = {q 6 Q \\q — pi\\ < Ri} and the half 
planes Sij are 



{q € R | 2q ■ (pi - p 3 ) > fa + Pj ). (pi - Pj )}. 



Provided Ri is twice as large as the maximum distance be- 
tween pi and the vertexes of W{pi,R{), all Voronoi neigh- 
bors of pi are within distance Ri from pi and the equality 
Vi = W(pi,Ri) holds. The minimum adequate sensing ra- 
dius is therefore Ri tmnl = 2max gg vi/(p ii i?, i min ) \\pi — q\\. We 
are now ready to state the following algorithm. 



A similar algorithm can be designed for a robotic agent 
with communication capabilities. The specifications go as 
in the previous algorithm, except for the fact that steps 2: 
and 7: are substituted by 

send ("request to reply", Pi{ti)) within radius Ri 
receive ("response", p^ from all agents within radius Ri 

Further, we have to require each agent to perform the fol- 
lowing event-driven task: if the ith agent receives at any 
time ti a "request to reply" message from the jth agent 
located at position pj , it executes 

send ("response", ^ (ti)) within radius \\pi(t) — Pj\\ 

We call this algorithm Adjust communication radius 

ALGORITHM. 

Next, we present an algorithm whose objective is to 
maintain the information about the Voronoi cell of the ith 
agent, and detect the presence of certain events. We con- 
sider only robotic agents with sensing capabilities. We call 
an agent active if it is moving and we assume that the ith 
agent can determine if any agent within radius Ri is active 
or not. Two events are of interest: (i) a Voronoi neighbor of 



8 



SUBMITTED AS A REGULAR PAPER TO IEEE TRANSACTIONS ON ROBOTICS AND AUTOMATION 



Name: Coverage behavior algorithm I 
Goal: distributed optimal agent location 

Requires: (i) Voronoi cell computation 

(ii) centroid and mass computation 

(iii) positive real So 

(iv) Adjust communication radius 

ALGORITHM 

For i £ {1, . . . , n}, ith agent performs at U = T^o = 0: 

0: P i (T i , ):=(p\(T i , ),...,pi l (T it0 )) 

0: compute Voronoi region Vi(Ti^) 

0: set Vi = Vi(T ifi ) and i? l = 2ma,x qeVi \\p t - q\\ 

For i G {1, ...,n}, the ith agent performs at time 
ti = Ti,k cither one of the following threads or both. 
For some Bi G N, we require that after Bi steps of 
the scheduling sequence, each of the threads has been 
executed at least once. 

[Information thread] 

i: run Adjust communication radius algorithm 
[Control thread] 

1: compute centroid Cy 4 and mass My i of Vi 
2: apply control pair {So, M Vi (C Vi - Pi{T i;k ))) 



the ith agent becomes active and (ii) a new active agent be- 
comes a Voronoi neighbor of the ith agent. In both cases, 
we require a trigger message "request recomputation" to 
an appropriate control algorithm that we shall present in 
the next section. Before presenting the algorithm, let us 
introduce the map weight that assigns to the state vector 
pi g ]g>ivxn a t U pi e (u)i f ...,w n ) G N" according to 



{3 if j G N{i) and j is active 
1 if j G M{i) and j is not active 




The algorithm is designed to run for times U G [to, to + Si\. 





Name: Monitoring algorithm 




Goal: Cell maintenance & event detection 




Requires: (i) sensor with radius Ri 




(ii) positive reals t , St 




(iii) Adjust sensing radius algo- 




rithm 


Local agent i performs for ti G [to, t +5t]: 


1 


initialize P l {to) and Vi(to), set w = weight (P l {to)) 


2 


while ti < to + St do 


3 


run Adjust sensing radius algorithm 


4 


if weighty (P i (ti)) >Wj+2 then 


5 


send ( "request recomputation" ) 


6 


set w = weight (P l (ti)) 


7 


end if 


8 


end while 



As a consequence of Theorem 3.1 and Corollary 3.1 
in pep , wc have the following result. 



C. Asynchronous distributed implementations of coverage 
control 



Let us now present two versions of Lloyd algorithm for 
the solution of the optimization problem ([j]) that can be 
implemented by an asynchronous distributed network of 
robotic agents. For simplicity, we assume that at time 
all clocks are synchronized (although they later can run at 
different speeds) and that each agent knows at the ex- 
act location of every other agent. The first algorithm is 
designed for robotic agents with communication capabil- 
ities, and requires the Adjust communication radius 
algorithm (while it does not require the Monitoring 
algorithm). 



Proposition IV. 1: Let Pq G Q n denote the initial sensors 
location. Let {Tj.} be the sequence in increased order of all 
the scheduling sequences of the agents of the network. As- 
sume mik{Tk — Tfc-i} > 0. Then, there exists a sufficiently 
small <5* > such that if < <5 < <5» , the Coverage be- 
havior algorithm I converges to the set of critical points 
of TLv, that is, the set of centroidal Voronoi configurations. 



Next, we focus on distributed asynchronous implementa- 
tions of Lloyd algorithm that take advantage of the special 
structure of the coverage problem. The following algorithm 
is designed for robotic agents with sensing capabilities, it 
requires the Monitoring and the Adjust sensing radius algo- 
rithms. Two advantages of this algorithm over the previous 
one are that there is no need for each agent to exactly go 
toward the centroid of its Voronoi cell nor to take a small 
step at each stage. 
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Name: Coverage behavior algorithm II 
Goal: distributed optimal agent location 

Requires: (i) Voronoi cell computation 

(ii) centroid computation 
(iii) Monitoring algorithm 



For i € {1, . . . , n}, ith agent performs at ij = T^o = 0: 



P^o) :=(piN,.,p;(r,,o)) 



compute Voronoi region T^(T^o) 
set Vi = Vi(T it0 ) and R t = 2max, 6 y i \\pi - q\\ 

For i G {1, . . . , n}, ith agent performs at ti = T^: 



1: choose < Sti < ti >m i n 

2: set s = Ti^, compute centroid Cy^s) 

3: choose Ui, with ui ■ (Cy i — Pi(s)) > 0, with strict 

inequality if Pi{s) ^ Cy i 
4: while ti < Ti : k + Sti do 

run Monitoring algorithm for (T^,^) 
while no warning do 

Pi = Ui 

end while 

set s = ti, compute centroid Cv i ( s ) 
10: choose Ui, with Ui • (Cy 4 — Pi{s)) > 0, with strict 

inequality if Pi(s) ^ Cv i 
11: end while 



Remark IV. 2: The control law Ui in step 7: can be de- 
fined via a saturation function. For instance, SR : R N — > 

x if ||x|| < 1 
x/\\x\\ if||x||>l 



SR(x) 



Then set Ui = SR(Cy i — Pi). 



Resorting to the discussion in Section III-B on the con 



vcrgcncc of the discrete Lloyd algorithms, one can prove 
that the Coverage behavior algorithm II veri fies p roperties 
(a) and (b). As a consequence of Proposition III.3| , we then 
have the following result. 

Proposition IV. 3: Let Pq G Q n denote the initial sensors 
location. The Coverage behavior algorithm II con- 
verges to the set of critical points of TCv, that is, the set of 
centroidal Voronoi configurations. 

V. Extensions and applications 

In this section we investigate various extensions and ap- 
plications of the algorithms proposed in the previous sec- 
tions. We extend the treatment to vehicles with passive 
dynamics and we also consider discrete-time implementa- 
tions of the algorithms for vehicles endowed with a local 
motion planner. Finally, we describe interesting ways of 
designing density functions to solve problems apparently 
unrelated to coverage. 

A. Variations on vehicle dynamics 

Here, we consider vehicles systems described by more 
general linear and nonlinear dynamical models. 

Coordination of vehicles with passive dynamics. We start 
by considering the extension of the control design to non- 
linear control systems whose dynamics is passive. Relevant 



examples include networks of vehicles and robots with gen- 
eral Lagrangian dynamics, as well as spatially invariant 
passive linear systems. Specifically, assume that for each 
i G {l,...,n}, the ith vehicle state includes the spatial 
variable pi, and that the vehicle's dynamics is passive with 
input Ui, output pi and storage function Si : Q — > M + . 
Furthermore, assume that the input preserving the zero 
dynamics manifold {pi = 0} is = 0. 

For such systems, we devise a proportional derivative 
(PD) control via, 



-k pIO pMv z (pi - CVi) - fcderivPi 



(12) 



where fc pr0 p and fcd er iv are scalar positive gains. The closed- 
loop system induced by this control law can be analyzed 
with the Lyapunov function 



1 



£ — T^propT^V + ^ ^ Si, 



i=l 

yielding the following result. 

Proposition V.l: For passive systems, the control 
law achieves asymptotic convergence of the sensors 
location to the set of centroidal Voronoi configurations. If 
this set is finite, then the sensors location converges to a 
centroidal Voronoi configuration. 

Proof: Consider the evolution of the function £, 



-£ — „ fcprop , 7~(-V + Si 



dt 2 plup dt 

1=1 

71 

< k P TopM Vt (pi~C Vi ) + PiUi = -fcdcriv ^ Pj < ■ 

i=l 

By LaSalle's principle, the sensors location converges to the 
largest invariant set contained in {pi = 0}. Given the as- 
sumption on the zero dynamics, we conclude that pi = Cy i 
for i G {1, . . . , n}, i.e., the largest invariant set corresponds 
to the set of centroidal Voronoi configurations. If this set 
is finite, LaSalle's principle also guarantees convergence to 
a specific centroidal Voronoi configuration. ■ 
In Figure || we illustrate the performance of the control 
law ( |l2| ) for vehicles with second-order dynamics pi = u^ . 




Fig. 3 

Coverage control for 32 vehicles with second order 
dynamics. The environment and Gaussian density function 
are as in Figure ^. The control gains are fcp RO p = 6 and 

^DERIV — 1- 



Coordination of vehicles with local controllers. Next, 
consider the setting where each vehicle has an arbitrary 



10 



SUBMITTED AS A REGULAR PAPER TO IEEE TRANSACTIONS ON ROBOTICS AND AUTOMATION 



dynamics and is endowed with a local feedback and feed- 
forward controller. The controller is capable of strictly de- 
creasing the distance to any specified position in Q in a 
specified period of time S. 

Assume the dynamics of the ith vehicle is described by 
Xi = fi(t,Xi,u), where Xi G R m denotes its state, and 
7T,; : M. mi — > Q is such that ni(xi) = p*. Assume also that 
for any p t ar g ot 6 Q and any x G R m \ 7rr 1 (p taxget ), there 
exists u(t, x(t) , parget) such that the solution x~i(t) of 

% = fi(t,Xi(t),u(t,Xi(t),p ta r g et)) , Xi(0) = Xq , 
Verifies \\iri(x t (t +6)) - Ptarget || < |Ki(Zi(*o)) ~ Ptarget ]|- 

Proposition V.2: Consider the following coordination al- 
gorithm. At time tk = k5, k G N, each vehicle computes 
Vi(tk) and Cvi(tk); then, for time t G [tk,tk+i{, the ve- 
hicle executes u(t, x(t), Cy i (tk))- For this closed-loop sys- 
tem, the sensors location converges to the set of centroidal 
Voronoi configurations. If this set is finite, then the sensors 
location converges to a centroidal Voronoi configuration. 
The proo f of this result readily follows from Proposi- 
tion Ill.j , sin ce the algorithm verifies properties (a) and 
(b) of Section |ibB . 

As an example, we consider a classic model of mobile 
wheeled dynamics, the unicycle model. Assume the ith 
vehicle has configuration (0i,Xi,yi) G SE(2) evolving ac- 
cording to 



Vi cos 6; 



where (w,-, v{) are the control inputs for vehicle i. Note that 
the definition of (9i,Vi) is unique up to the discrete action 
(9i, Vi) i ► (9i + 7r, —Vi). Given a target point ptarget, we use 
this symmetry to require the equality (cos#i, sin#i) • fa — 
Ptarget) < for all time t. Should the equality be violated 
at some time t = to, we shall redefine 9i(t^) = 9i(t^) + n 
and Vi as — Vi from time t — to onwards. 

Following the approach in |52j , consider the control law 



_ „, (-sinflj, C0S6>^) ■ fa - Ptarget) 

uii — ZK pr0 p arctan — ; — — — — 

(COS Vi , Sin t ) ■ fa— ^target ) 
Vi = -fc pr op( cos ^: Sin^) • fa - Ptarget), 



where fc prop is a positive gain. This feedback law differs 
from the original stabilizing strategy in J5^] only in the fact 
that no final angular position is preferred. One can prove 
that pi — (xi, yi) is guaranteed to monotonically approach 
the target position ptarget when run over an infinite time 
horizon. We illustrate the performance of the proposed 
algorithm in Figure ^. 

B. Geometric patterns and formation control 

Here we suggest the use of decentralized coverage algo- 
rithms as formation control algorithms, and we present 
various density functions that lead the multi-vehicle net- 
work to predetermined geometric patterns. In particular, 
we present simple density functions that lead to segments, 
ellipses, polygons, or uniform distributions inside convex 
environments. 



Consider a planar environment, let k be a large positive 
gain, and denote q = (x, y) G Q C K 2 . Let a,b,c be real 
numbers, consider the line ax + by + c = 0, and define the 
density function 

<hne(q) = cxp(-k(ax + by + c) 2 ). 

Similarly, let (x c , y c ) be a reference point in M 2 , let a, b, r be 
positive scalars, consider the ellipse a(x — x c ) 2 +b(y — y c ) 2 = 
r 2 , and define the density function 

<?Wipsc(<7) = exp ( - k(a(x - x c ) 2 + b(y - y c ) 2 - r 2 ) 2 ). 

We illustrate this density function in Figure |^. During 
the simulations, we observed that the convergence to the 
desired pattern was rather slow. 




Fig. 5 

Coverage control for 32 vehicles with ^ellipse- The 

PARAMETER VALUES ARE: k = 500, a = 1.4, b = .6, x c = y c = 0, 
r 2 = .3, AND fcp R op = 1. 



Finally, define the smooth ramp function SR^x) = 
a;(arctan(fo;)/7r + (1/2)), and the density function 

0disk(<?) = cxp(-fc SR<>(a(.T - x c ) 2 + b(y - y c ) 2 - r 2 )). 

This density function leads the multi-vehicle network to 
obtain a uniform distribution inside the ellipsoidal disk 



a(x - x c ) 2 + b{y - y c ) 2 
function in Figure O. 



< r 1 



We illustrate this density 




Fig. 6 

Coverage control for 32 vehicles to an ellipsoidal disk. The 
density function parameters are the same as in figure ^, and 
£ — 10, fcpRop — 1- 



It appears straightforward to generalize these types of 
density functions to the setting of arbitrary curves or 
shapes. The proposed algorithms are to be contrasted with 
the classic approach to formation control based on rigidly 
encoding the desired geometric pattern. One disadvantage 
of the proposed approach is the requirement for a careful 
numerical computation of Voronoi diagrams and centroids. 
We refer to [jl6| for previous work on algorithms for ge- 
ometric patterns, and to 0, E3] for formation control 
algorithms. 
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| / > 



VTi 



Fig. 4 

Coverage control for 16 vehicles with mobile wheeled dynamics. The environment and Gaussian density function are as in 

Figure bl and fc PB op = 3. 



VI. Conclusions 

We have presented a novel approach to coordination al- 
gorithms for multi- vehicle networks. The scheme can be 
thought of as an interaction law between agents and as 
such it is implcmcntable in a distributed asynchronous 
fashion. Numerous extensions appear worth pursuing. We 
plan to investigate the setting of non-convex environments 
and non-isotropic sensors. We are currently implementing 
these algorithms on a network of all-terrain vehicles. Fur- 
thermore, we plan to extend the algorithms to provide col- 
lision avoidance guarantees and to vehicle dynamics which 
are not locally controllable. 

Acknowledgments 

This work was supported by NSF Grant CMS-0100162, 
ARO Grant DAAD 190110716, and DARPA / AFOSR 
MURI Award F49620-02-1-0325. 

References 

[1] C R. Wcisbin, J. Blitch, D. Lavcry, E. Krotkov, C Shoemaker, 
L. Matthies, and G. Rodriguez, "Miniature robots for space and 
military missions," IEEE Robotics & Automation Magazine, vol. 
6, no. 3, pp. 9-18, 1999. 

[2] E. Krotkov and J. Blitch, "The Defense Advanced Research 
Projects Agency (DARPA) tactical mobile robotics program," 
International Journal of Robotics Research, vol. 18, no. 7, pp. 
769-76, 1999. 

[3] P. E. Rybski, N. P. Papanikolopoulos, S. A. Stoeter, D. G. 
Krantz, K. B. Yesin, M. Gini, R. Voyles, D. F. Hougen, B. Nel- 
son, and M. D. Erickson, "Enlisting rangers and scouts for re- 
connaissance and surveillance," IEEE Robotics & Automation 
Magazine, vol. 7, no. 4, pp. 14-24, 2000. 

[4] T. B. Curtin, J. G. Bellingham, J. Catipovic, and D. Webb, 
"Autonomous oceanographic sampling networks," Oceanogra- 
phy, vol. 6, no. 3, pp. 86-94, 1993. 

[5] R. M. Turner and E. H. Turner, "Organization and reorganiza- 
tion of autonomous oceanographic sampling networks," in IEEE 
Int. Conf. on Robotics and Automation, Leuven, Belgium, May 
1998, pp. 2060-7. 

[6] E. Eberbach and S. Phoha, "SAMON: communication, coop- 
eration and learning of mobile autonomous robotic agents," in 
Proceedings 11th International Conf. on Tools with Artificial 
Intelligence (TAI), Chicago, IL, Nov. 1999, pp. 229-36. 

[7] A. Okabe, B. Boots, and K. Sugihara, "Nearest neighbourhood 
operations with generalized Voronoi diagrams: a review," In- 
ternational Journal of Geographical Information Systems, vol. 
8, no. 1, pp. 43-71, 1994. 



[8] Z. Drezner, Ed., Facility Location: A Survey of Applications 
and Methods, Springer Series in Operations Research. Springer 
Verlag, New York, NY, 1995. 

[9] A, Suzuki and A. Okabe, "Using Voronoi diagrams," In Drezner 
[§, pp. 103-118. 

[101 Okabe and A. Suzuki, "Locational optimization problems 
solved through Voronoi diagrams," European Journal of Oper- 
ational Research, vol. 98, no. 3, pp. 445-56, 1997. 
[Ill A. Okabe, B. Boots, K. Sugihara, and S. N. Chiu, Spatial Tessel- 
lations: Concepts and Applications of Voronoi Diagrams, Wiley 
Series in Probability and Statistics. John Wiley & Sons, New 
York, NY, second edition, 2000. 
[12] Q. Du, V. Faber, and M. Gunzburger, "Centroidal Voronoi 
tessellations: applications and algorithms," SIAM Review, vol. 
41, no. 4, pp. 637-676, 1999. 
[131 S. P. Lloyd, "Least squares quantization in PCM," IEEE Trans- 
actions on Information Theory, vol. 28, no. 2, pp. 129—137, 1982, 
Presented as Bell Laboratory Technical Memorandum at a 1957 
Institute for Mathematical Statistics meeting. 
[14] R. M. Gray and D. L. Neuhoff, "Quantization," IEEE Trans- 
actions on Information Theory, vol. 44, no. 6, pp. 2325-2383, 
1998, Commemorative Issue 1948-1998. 
[15] H. Yamaguchi and T. Arai, "Distributed and autonomous 
control method for generating shape of multiple mobile robot 
group," in IEEE/RSJ Int. Conf. on Intelligent Robots & Sys- 
tems, Munich, Germany, Sept. 1994, pp. 800-807. 
[16] K. Sugihara and I. Suzuki, "Distributed algorithms for formation 
of geometric patterns with many mobile robots," Journal of 
Robotic Systems, vol. 13, no. 3, pp. 127—39, 1996. 
[17] T. Balch and R. Arkin, "Behavior-based formation control for 
multirobot systems," IEEE Transactions on Robotics and Au- 
tomation, vol. 14, no. 6, pp. 926—39, 1998. 
[181 M. Egerstedt and X. Hu, "Formation constrained multi-agent 
control," IEEE Transactions on Robotics and Automation, vol. 
17, no. 6, pp. 947-51, 2001. 
[19] J. P. Dcsai, J. P. Ostrowski, and V. Kumar, "Modeling and 
control of formations of nonholonomic mobile robots," IEEE 
Transactions on Robotics and Automation, vol. 17, no. 6, pp. 
905-8, 2001. 

[20] P. Tabuada, G. Pappas, and P. Lima, "Feasible formations of 

multi-agent systems," Automatica, 2002, Submitted. 
[211 R- Olfati-Saber and R. M. Murray, "Graph rigidity and dis- 
tributed formation stabilization of multi-vehicle systems," in 
IEEE Conf. on Decision and Control, Las Vegas, NV, 2002, To 
appear. 

[22] C. Tomlin, G. J. Pappas, and S. S. Sastry, "Conflict resolution 
for air traffic management: a study in multiagent hybrid sys- 
tems," IEEE Transactions on Automatic Control, vol. 43, no. 
4, pp. 509-21, 1998. 

[23] E. Frazzoli, Z. H. Mao, J. H. Oh, and E. Feron, "Aircraft conflict 
resolution via semi-definite programming," AIAA Journal of 
Guidance, Control, and Dynamics, vol. 24, no. 1, pp. 79-86, 
2001. 



SUBMITTED AS A REGULAR PAPER TO IEEE TRANSACTIONS ON ROBOTICS AND AUTOMATION 



H. Choset, "Coverage for robotics - a survey of recent results," 
Annals of Mathematics and Artificial Intelligence, vol. 31, pp. 
113-126, 2001. 

R. Bachmayer and N. Ehrich Leonard, "Vehicle networks for 
gradient descent in a sampled environment," in IEEE Conf. on 
Decision and Control, Dec. 2002, To appear. 
A. Jadbabaie, J. Lin, and A. S. Morse, "Coordination of groups 
of mobile autonomous agents using nearest neighbor rules," 
IEEE Transactions on Automatic Control, July 2002, To ap- 
pear. 

E. Klavins, "Communication complexity of multi-robot sys- 
tems," in Workshop on Algorithmic Foundations of Robotics, 
Nice, France, Dec. 2002, Submitted. 

J. Cortes, S. Martinez, T. Karatas, and F. Bullo, "Coverage 
control for mobile sensing networks," in IEEE Int. Conf. on 
Robotics and Automation, Arlington, VA, May 2002, pp. 1327- 
1332. 

J. Cortes, S. Martinez, T. Karatas, and F. Bullo, "Coverage 
control for mobile sensing networks: variations on a theme," in 
Mediterranean Conference on Control and Automation, Lisbon, 
Portugal, July 2002, Electronic proceedings. 
R. A. Brooks, "A robust layered control-system for a mobile 
robot," IEEE Journal of Robotics and Automation, vol. 2, no. 

I, pp. 14-23, 1986. 

C. W. Reynolds, "Flocks, herds, and schools: A distributed 
behavioral model," Computer Graphics, vol. 21, no. 4, pp. 25- 
34, 1987. 

R. C. Arkin, Behavior-Based Robotics, Cambridge University 
Press, New York, NY, 1998. 

M. S. Fontan and M. J. Mataric, "Territorial multi-robot task 
division," IEEE Transactions on Robotics and Automation, vol. 
14, no. 5, pp. 815-822, 1998. 

A. C. Schultz and L. E. Parker, Eds., Multi-Robot Systems: 
From Swarms to Intelligent Automata, Washington, DC, June 
2002. Kluwer Academic Publishers, Proceedings from the 2002 
NRL Workshop on Multi-Robot Systems. 

T. Balch and L. E. Parker, Eds., Robot Teams: From Diversity 
to Polymorphism, A. K. Peters Ltd., 2002. 

L. E. Parker, "Distributed algorithms for multi-robot observa- 
tion of multiple moving targets," Autonomous Robots, vol. 12, 
no. 3, pp. 231-55, 2002. 

A. Howard, Maja J. Mataric, and G. S. Sukhatmc, "Mobile 
sensor network deployment using potential fields: A distributed 
scalable solution to the area coverage problem," in Proceedings 
of the 6th International Conference on Distributed Autonomous 
Robotic Systems (DARS02), Fukuoka, Japan, 2002, pp. 299-308. 
N. A. Lynch, Distributed Algorithms, Morgan Kaufmann Pub- 
lishers, San Mateo, CA, 1997. 

G. Tel, Introduction to Distributed Algorithms, Cambridge Uni- 
versity Press, New York, NY, second edition, 2001. 
J. N. Tsitsiklis, D. P. Bertsekas, and M. Athans, "Distributed 
asynchronous deterministic and stochastic gradient optimization 
algorithms," IEEE Transactions on Automatic Control, vol. 31, 
no. 9, pp. 803-12, 1986. 

D. P. Bertsekas and J. N. Tsitsiklis, Parallel and Distributed 
Computation: Numerical Methods, Athena Scientific, 1997. 

S. H. Low and D. E. Lapscy, "Optimization flow control I: Ba- 
sic algorithm and convergence," IEEE/ACM Transactions on 
Networking, vol. 7, no. 6, pp. 861-74, 1999. 

G. Leibon and D. Lctscher, "Delaunay triangulations and 
Voronoi diagrams for Riemannian manifolds," in Proceedings 
of the Sixteenth Annual Symposium on Computational Geome- 
try (Hong Kong, 2000), New York, 2000, pp. 341-349, ACM. 
R. Klein, Concrete and abstract Voronoi diagrams, vol. 400 of 
Lecture Notes in Computer Science, Springer Verlag, New York, 
NY, 1989. 

A. Suzuki and Z. Drezner, "The p-center location problem in an 

area," Location Science, vol. 4, no. 1/2, pp. 69—82, 1996. 

J. Cortes and F. Bullo, "Distributed Lloyd flows for disk covering 

and sphere packing problems," 2002, Preprint. 

C. Cattani and A. Paoluzzi, "Boundary integration over linear 

polyhedra," Computer- Aided Design, vol. 22, no. 2, pp. 130-5, 

1990. 

J. Gao, L. J. Guibas, J. Hershbcrger, Li Zhang, and An Zhu, 
"Geometric spanner for routing in mobile networks," in ACM 
International Symposium on Mobile Ad-hoc Networking & Com- 
puting, Long Beach, CA, Oct. 2001, pp. 45-55. 
X.-Y. Li and P.-J. Wan, "Constructing minimum energy mobile 



wireless networks," ACM Journal of Mobile Computing and 
Communication Survey, vol. 5, no. 4, 2001. 
[50] S. Meguerdichian, S. Slijepcevic, V. Karayan, and M. Potkinjak, 
"Localized algorithms in wireless ad- hoc networks: Location dis- 
covery and sensor exposure," in ACM International Symposium 
on Mobile Ad-hoc Networking & Computing, Long Beach, CA, 
Oct. 2001. 

[51] M. Cao and C. Hadjicostis, "Distributed algorithms for Voronoi 
diagrams and application in ad-hoc networks," Preprint, Oct. 
2002. 

[52] A. Astolfi, "Exponential stabilization of a wheeled mobile robot 
via discontinuous control," ASME Journal on Dynamic Sys- 
tems, Measurement, and Control, vol. 121, no. 1, pp. 121-7, 
1999. 

[53] H. K. Khalil, Nonlinear Systems, Prentice Hall, Englcwood 

Cliffs, NJ, second edition, 1995. 
[54] D. G. Luenbcrger, Linear and Nonlinear Programming, 

Addison- Wesley, Reading, Massachusetts, second edition, 1984. 

VII. Appendix 

In this section we collect some relevant facts on descent 
flows both in the continuous and in the discrete-time set- 
tings. We do this following |53| and pM, respectively. We 



include Proposition VII. 4 as we are unable to locate it in 



the linear and nonlinear programming literature. 

Continuous-time descent flows 

Consider the differential equation x = X(x), where X : 
D C M. N — ► is locally Lipschitz and D is an open 
connected set. A set M is said to be (positively) invariant 
with respect to X if x(0) £ M implies x{t) £ M, for all 
t £ R (resp. t > 0). A descent function for X on f2, 
V : D — ► R, is a continuously differentiable function such 
that CxV < on f2. We denote by E the set of points in 
f2 where CxV(x) = and by M be the largest invariant 
set contained in E. Finally, the distance from a point x to 
a set M is defined as dist(x, M) = inf pg jvf \\ x — T>\\- 

Lemma VII. 1 (LaSalle's principle) Let fl C D be a com- 
pact set that it is positively invariant with respect to X. 
Let x(0) £ M and be an accumulation point of x(t). 
Then x* £ M and dist(x(i), M) -> as t -> oo. 
The following corollary is Exercise 3.22 in fl53| . 

Corollary VII. 2: If the set M is a finite collection of 
points, then the limit of x(t) exists and equals one of them. 

Discrete-time descent flows 

Let X be a subset of M. N . An algorithm T is a continuous 
mapping from X to X. A set C is said to be positively 
invariant with respect to T if xo £ C implies T{xq) £ C . 
A point x* is said to be a fixed point of T if T(x*) = X*. 
We denote the set of fixed points of T by T. A descent 
function for TonC, Z : X — > R+, is any nonnegative real- 
valued continuous function satisfying Z(T{x)) < Z(x) for 
x £ C, where the inequality is strict if x L. Typically, 
Z is the objective function to be minimized, and T reflects 
this goal by yielding a point that reduces (or at least does 
not increase) Z. 

Lemma VII. 3 (Global convergence theorem) Let C C 
X be a compact set that it is positively invariant with re- 
spect to T. Let xq £ C and denote x rn = T(x m -i), m > 1. 
Let x* be an accumulation point of the sequence {x m } m >i. 
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Then x* £ T, dist(x m ,r) — ► and Z(x m ) — ► Z(x*) as 
m — > oo. 

Proposition VII. 4: If the set T is a finite collection of 
points, then {x m } converges and equals one of them. 

Proof: Let x* be an accumulation point of {x m } 
and assume the whole sequence does not converge to it. 
Then, there exists an e > such that for all too, there 
is a ml > too such that \\x m ' ~ #*|| > e. Let d be the 
minimum of all the distances between the points in L. Fix 
e' = min{d/2, e}. Since T is continuous and F is finite, 
there exists S > such that ||x — z\\ < S, with z £ F, 
implies ||T(x) — z|| < e' (that is, for each z £ T, there 
exists such S(z), and we take the minimum over T). 

Now, since dist(x m , T) — > 0, there exists toi such that for 
all to > Toi, dist(x m , T) < S. Also, we know that there is a 
subsequence of {x m } which converges to x*, let us denote 
it by {x mk }k>i- For 5, there exists such that for all 
k > k\, we have \\x mk — x*|| < 5. 

Let Too = maxjmi, to/^ }. Take fc such that to^ > mo 
Then, 

\\x mk+ i - x*|| = ||T(x m J - x*|| < e' . (13) 

Now we are going to prove that ||x rrafc +i — x*|| < 5. If 
d/2 < S, then this claim is straightforward, since e' < d/2. 
If d/2 > S, suppose that ||x 7Tifc +i — x* || > 5. Since to^ + 1 > 
fno > fni, then dist(x„ ifc +i, T) < 5. Therefore, there exists 
y £ r such that ||x mfc +i — y\\ < S. Necessarily, y ^ x*. 
Now, by the triangle inequality, ||x — y\\ < \\x — x mk+ i\\ + 
\\x mk +i - y\\. Then, 

\\x„ lk+1 -x*|| > ||x-y|| - ||x mfc+1 -y|| > d- 6 > d/2, 

which contradicts (|lj). Therefore, ||x TO( .+i — x*|| < 5. This 
argument can be iterated to prove that for all to > too, wc 
have ||x m — x*|| < S. Let us take now ml > too such that 
||x m / — x*|| > e. Since to' — 1 > too, we have ||x m '_i — x» || < 
(5, and therefore 

\\x m > - x*\\ = ||T(x m /_i) - x*|| < e < e , 

which is a contradiction. Therefore, {x TO } converges to x*. 



